#R version 4.4.3 (2025-02-28 ucrt)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 11 x64 

library(dplyr) #1.1.4
library(ggplot2) #3.5.1
library(readr) #2.1.5
library(haven) #2.5.4
library(estimatr) #1.0.6
library(broom.mixed) #0.2.9.6
library(lme4) #1.1-37
library(countrycode) #1.6.1
library(forcats) #1.0.0
library(stargazer) #5.2.3
library(tibble) #3.2.1
library(modelsummary) #2.3.0
library(clubSandwich) #0.6.0

`%notin%` <- function(x,y) !(x %in% y) 

setwd("C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/cps/replication")

df <- read_csv("data/gwf_panel.csv") |>  
  select(countryname, year, v2xps_party, v2x_polyarchy, gwf_duration, gwf_regime, gwf_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm) |>  
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown"), 
         reactive = case_when(ruling == 1 ~ 0, ruling == 0 ~ 1)) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.nan(.) | gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.na(.), 0, .)) |>  
  filter(gwf_prior %in% c("Party", "Military", "Personal") & gwf_regime == "democracy") |>  
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) |>  
  mutate(gwf_prior = factor(gwf_prior, levels = c("Party", "Personal", "Military"))) |>  
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) |> 
  group_by(ccode) |> 
  mutate(ldv = lag(v2xps_party), reactive = case_when(ruling == 1 ~ 0, ruling == 0 ~ 1))

#Figure 1----
gwf <- df
mp <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = gwf, se_type = "CR0", clusters = ccode) |>
  tidy() |>
  mutate(Sample = "Pooled") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

df1 <- filter(gwf, elections == "None")

m0 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "CR0", clusters = ccode) |>
  tidy() |>
  mutate(Sample = "None") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

df1 <- filter(gwf, elections == "First")

m1 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "CR0", clusters = ccode) |>
  tidy() |>
  mutate(Sample = "First") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

df1 <- filter(gwf, elections == "Second")

m2 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "CR0", clusters = ccode) |>
  tidy() |>
  mutate(Sample = "Second") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

df1 <- filter(gwf, elections == "Third-plus")

m3 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "CR0", clusters = ccode) |>
  tidy() |>
  mutate(Sample = "Third-plus") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

#fig1 left hand side
rbind(m0, m1, m2, m3, mp) |>
  mutate(Sample = factor(Sample, levels = c("None", "First", "Second", "Third-plus", "Pooled"))) |>
  mutate(term = recode(term, "pi1" = "Incumbent PI", "prior_pi" = "Prior PI", "asp_dummy" = "ASP", 
                       "gwf_priorParty" = "Party", "gwf_priorPersonal" = "Personal", "gwf_priorMilitary" = "Military")) |>
  filter(term %in% c("Incumbent PI", "Prior PI", "ASP")) |>
  ggplot(aes(x = Sample, y = estimate, ymin = conf.low, ymax = conf.high, color = term, shape = term)) +
  geom_pointrange(position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(y = estimate - 100), linetype = "solid", color = NA) +  # Setting color to NA to remove legend
  geom_point(aes(y = estimate - 100), shape = NA, color = NA) +  # Setting color and shape to NA to remove legend
  theme_bw() +
  xlab("") +
  ylab("") +
  ylim(-0.35, 0.65) +
  scale_color_grey(name = "Variable", start = 0, end = 0.8) +
  scale_shape_manual("Variable", values = c(15, 16, 17)) + 
  theme(legend.position = "bottom")  +
  ggtitle("OLS Models") +
  theme(plot.title = element_text(hjust = 0.5))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_d_1_1.jpg", dpi = 500, width = 5, height = 5)

#Table 2 CRSE----
gwf2 <- gwf |> select(v2xps_party, v2x_polyarchy, lgdp, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()
gwf3 <- gwf |> select(countryname, year, v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, ccode) |> na.omit()
gwf4 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, gwf_pduration, ccode) |> na.omit()
gwf4 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, gwf_pduration, regime_elecs, ccode) |> na.omit()
gwf5 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, gwf_number, fsu, ussr_sat, av_dm, ccode) |> na.omit()
gwf6 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, gwf_number, fsu, ussr_sat, av_dm, ccode) |> na.omit()
gwf7 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()
gwf8 <- gwf |> select(ldv, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()
gwf9 <- gwf |> select(year, ccode, ldv, gwf_prior, pi1, prior_pi, reactive, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()

m1 <- lmer(v2xps_party ~ + (1 | ccode), data = gwf, REML=FALSE) #

m2 <- lmer(v2xps_party ~ v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf2, REML=FALSE) #

m3 <- lmer(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + (1 | ccode), data = gwf3, REML=FALSE) # 

m4 <- lmer(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + gwf_pduration + (1 | ccode), data = gwf4, REML=FALSE) #

m5 <- lmer(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + av_dm + (1 | ccode), data = gwf5, REML=FALSE)

m7 <- lmer(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf7, REML=FALSE)

m8 <- lm(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = gwf7) 
sd(residuals(m8))

m9 <- lm(ldv ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = gwf8) 
sd(residuals(m9))

m10 <- lm(ldv ~ gwf_prior + pi1 + prior_pi + reactive + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + regime_elecs + av_dm, data = gwf9) 
sd(residuals(m10))

modellist <- list("ML1" = m1, #just baseline
                  "ML2" = m2, #shows basline with time-variant indicators
                  "ML3" = m3, #just time invariant - prior duration 
                  "ML4" = m4, #add in prior duration, not much change 
                  "ML5" = m5, #full minus prior duration & # of elections, prior_pi drops 
                  "ML7" = m7,
                  "LM" = m8,
                  "LDV" = m9,
                  "LDV" = m10)

crse1 <- clubSandwich::vcovCR(m1, cluster = gwf$ccode, type = "CR0")
crse2 <- clubSandwich::vcovCR(m2, cluster = gwf2$ccode, type = "CR0")
crse3 <- clubSandwich::vcovCR(m3, cluster = gwf3$ccode, type = "CR0")
crse4 <- clubSandwich::vcovCR(m4, cluster = gwf4$ccode, type = "CR0")
crse5 <- clubSandwich::vcovCR(m5, cluster = gwf5$ccode, type = "CR0")
crse7 <- clubSandwich::vcovCR(m7, cluster = gwf7$ccode, type = "CR0")
crse8 <- clubSandwich::vcovCR(m8, cluster = gwf7$ccode, type = "CR0")
crse9 <- clubSandwich::vcovCR(m9, cluster = gwf8$ccode, type = "CR0")
crse10 <- clubSandwich::vcovCR(m10, cluster = gwf9$ccode, type = "CR0")

selist <- list(crse1, crse2, crse3, crse4, crse5, crse7, crse8, crse9, crse10)

modelsummary(modellist, stars = TRUE, gof_omit = "R2|AIC|RM|F|Log.Lik", vcov = selist, coef_map = c("pi1" = "Incumbent PI", "prior_pi" = "Prior PI", "asp_dummy" = "ASP",
                                                                                     "reactive" = "Reactive ASP",
                                                                                     "gwf_priorPersonal" = "Personalist", "gwf_priorMilitary" = "Military",
                                                                                     "gwf_pduration" = "Prior Duration", "regime_elecs" = "Elections", 
                                                                                     "SD (Intercept ccode)" = "SD: Country", "SD (Observations)" = "SD: Residuals"), output = "latex", fmt = 3)


#Figure 2 ----
df <- read_csv("data/gwf_panel.csv") |>  
  select(countryname, year, v2xps_party, v2x_polyarchy, gwf_duration, gwf_regime, gwf_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm) |>  
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown"), 
         reactive = case_when(ruling == 1 ~ 0, ruling == 0 ~ 1)) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.nan(.) | gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.na(.), 0, .)) |>  
  filter(gwf_prior %in% c("Party", "Military", "Personal") & gwf_regime == "democracy") |>  
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) |>  
  mutate(gwf_prior = factor(gwf_prior, levels = c("Party", "Personal", "Military"))) |>  
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus")))

# Function to run the model with a specific baseline year
pi_models <- function(data, baseline_year) {
  # Relevel gwf_duration to set baseline year
  data <- data |> mutate(gwf_duration = relevel(as.factor(gwf_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm_robust(v2xps_party ~ gwf_prior + pi1 * gwf_duration + prior_pi + asp_dummy, data = data, se_type = "CR0", clusters = ccode)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> 
    filter(term == "pi1")
  
  return(coefficients)
}
prior_models <- function(data, baseline_year) {
  # Relevel gwf_duration to set baseline year
  data <- data %>% mutate(gwf_duration = relevel(as.factor(gwf_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm_robust(v2xps_party ~ gwf_prior + prior_pi * gwf_duration + pi1 + asp_dummy, data = data, se_type = "CR0", clusters = ccode)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> 
    filter(term == "prior_pi")
  
  return(coefficients)
}
asp_models <- function(data, baseline_year) {
  # Relevel gwf_duration to set baseline year
  data <- data %>% mutate(gwf_duration = relevel(as.factor(gwf_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm_robust(v2xps_party ~ gwf_prior + asp_dummy * gwf_duration + pi1 + prior_pi, data = data, se_type = "CR0", clusters = ccode)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> 
    filter(term == "asp_dummy")
  
  return(coefficients)
}

# Create an empty list to store the results
pi_results <- list()
prior_results <- list()
asp_results <- list()

# Loop through the first 10 years to use each as the baseline
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- pi_models(df, baseline_year = as.numeric(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  pi_results[[as.numeric(year)]] <- coefficients
}
#prior models
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- prior_models(df, baseline_year = as.character(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  prior_results[[as.character(year)]] <- coefficients
}
#asp models
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- asp_models(df, baseline_year = as.character(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  asp_results[[as.character(year)]] <- coefficients
}

# Combine all results into a tidy dataframe
prior <- bind_rows(prior_results) 
pi <- bind_rows(pi_results) 
asp <- bind_rows(asp_results) 

bind_rows(prior, pi, asp) |> 
  mutate(term = case_when(term == "pi1" ~ "Incumbent PI", term == "prior_pi" ~ "Prior PI", term == "asp_dummy" ~ "ASP")) |> 
  ggplot(aes(x = as.factor(baseline_year), y = estimate, ymin = conf.low, ymax = conf.high, color = term, shape = term)) +
  geom_pointrange(show.legend = TRUE, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(y = estimate - 100), linetype = "solid", color = NA) +  # Setting color to NA to remove legend
  geom_point(aes(y = estimate - 100), shape = NA, color = NA) +  # Setting color and shape to NA to remove legend
  theme_bw() +
  xlab("") +
  ylab("") +
  ylim(-0.55, 0.55) +
  scale_color_grey(name = "Variable", start = 0, end = 0.8) +
  scale_shape_manual("Variable", values = c(15, 16, 17)) + 
  theme(legend.position = "bottom") +
  guides(color = guide_legend(order = 1), shape = guide_legend(order = 1))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_d_2_1.jpg", dpi = 500, width = 5, height = 5) 

# Function to run the model with a specific baseline year - without ASP------
pi_models <- function(data, baseline_year) {
  # Relevel gwf_duration to set baseline year
  data <- data |> mutate(gwf_duration = relevel(as.factor(gwf_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm_robust(v2xps_party ~ gwf_prior + pi1 * gwf_duration + prior_pi, data = data, se_type = "CR0", clusters = ccode)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> #filter(is.na(group)) |> 
    #mutate(ses = sqrt(diag(vcovCR(model, cluster = gwf$ccode, type = "CR2"))), model = "None", conf.low = estimate - (qnorm(0.975)*ses), conf.high = estimate + (qnorm(0.975)*ses)) |> 
    filter(term == "pi1")
  
  return(coefficients)
}
prior_models <- function(data, baseline_year) {
  # Relevel gwf_duration to set baseline year
  data <- data |> mutate(gwf_duration = relevel(as.factor(gwf_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm_robust(v2xps_party ~ gwf_prior + prior_pi * gwf_duration + pi1, data = data, se_type = "CR0", clusters = ccode)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> 
    filter(term == "prior_pi")
  
  return(coefficients)
}

# Create an empty list to store the results
pi_results <- list()
prior_results <- list()

# Loop through the first 10 years to use each as the baseline
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- pi_models(df, baseline_year = as.character(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  pi_results[[as.character(year)]] <- coefficients
}
#prior models
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- prior_models(df, baseline_year = as.character(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  prior_results[[as.character(year)]] <- coefficients
}

# Combine all results into a tidy dataframe
prior <- bind_rows(prior_results) 
pi <- bind_rows(pi_results) 

bind_rows(prior, pi) |> 
  mutate(term = case_when(term == "pi1" ~ "Incumbent PI", term == "prior_pi" ~ "Prior PI")) |> 
  ggplot(aes(x = as.factor(baseline_year), y = estimate, ymin = conf.low, ymax = conf.high, color = term, shape = term)) +
  geom_pointrange(show.legend = TRUE, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(y = estimate - 100), linetype = "solid", color = NA) +  # Setting color to NA to remove legend
  geom_point(aes(y = estimate - 100), shape = NA, color = NA) +  # Setting color and shape to NA to remove legend
  theme_bw() +
  xlab("") +
  ylab("") +
  ylim(-0.55, 0.55) +
  scale_color_grey(name = "Variable", start = 0.5, end = 0.8) +
  scale_shape_manual("Variable", values = c(16, 17)) + 
  theme(legend.position = "bottom") +
  guides(color = guide_legend(order = 1), shape = guide_legend(order = 1))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_d_2_2.jpg", dpi = 500, width = 5, height = 5) 

#reactive*pi1####
library(interplot) #0.2.3
library(marginaleffects) #0.25.1
lm1 <- lm(v2xps_party ~ gwf_prior + reactive*pi1 + prior_pi + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = df)
me1 <- plot_slopes(lm1, variables = "reactive", condition = "pi1", vcov = ~ccode) 
dme1 <- me1$plot_env$dat

interplot(m = lm1, var1 = "reactive", var2 = "pi1", hist = TRUE, adjCI = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of a Reactive ASP on \n Party Institutionalization") +
  xlab("Incumbent PI") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold")) +  
  ylim(-0.2, 0.5) +
  geom_ribbon(data = dme1, aes(x = pi1, y = estimate, ymin = conf.low, ymax = conf.high), alpha = 0.1) 

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_d_3_1.jpg", dpi = 500, width = 5, height = 5) 

lm2 <- lm(v2xps_party ~ gwf_prior + reactive*pi1 + prior_pi + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = df)
me2.1 <- plot_slopes(lm2, variables = "pi1", condition = "reactive") 
dme2.1 <- me2.1$plot_env$dat

me2.2 <- plot_slopes(lm1, variables = "pi1", condition = "reactive", vcov = ~ccode) 
dme2.2 <- me2.2$plot_env$dat

ggplot() +
  geom_pointrange(data = dme2.1, aes(x = reactive, y = estimate, ymin = conf.low, ymax = conf.high), linewidth = 1.2) +
  geom_pointrange(data = dme2.2, aes(x = reactive, y = estimate, ymin = conf.low, ymax = conf.high), linewidth = 0.5) +
  ylim(-0.25, 0.3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of Incumbent PI on \n Party Institutionalization") +
  xlab("Reactive ASP") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold"))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_d_3_2.jpg", dpi = 500, width = 5, height = 5) 

#reactive*prior_pi####
lm1 <- lm(v2xps_party ~ gwf_prior + reactive*prior_pi + pi1 + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = df)
me1 <- plot_slopes(lm1, variables = "reactive", condition = "prior_pi", vcov = ~ccode) 
dme1 <- me1$plot_env$dat

interplot(m = lm1, var1 = "reactive", var2 = "prior_pi", hist = TRUE, adjCI = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of a Reactive ASP on \n Party Institutionalization") +
  xlab("Prior PI") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold")) +  
  ylim(-0.35, 0.6) +
  geom_ribbon(data = dme1, aes(x = prior_pi, y = estimate, ymin = conf.low, ymax = conf.high), alpha = 0.1) 

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_d_4_1.jpg", dpi = 500, width = 5, height = 5) 

lm2 <- lm(v2xps_party ~ gwf_prior + reactive*prior_pi + pi1 + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = df)
me2.1 <- plot_slopes(lm2, variables = "prior_pi", condition = "reactive") 
dme2.1 <- me2.1$plot_env$dat

me2.2 <- plot_slopes(lm2, variables = "prior_pi", condition = "reactive", vcov = ~ccode) 
dme2.2 <- me2.2$plot_env$dat

ggplot() +
  geom_pointrange(data = dme2.1, aes(x = reactive, y = estimate, ymin = conf.low, ymax = conf.high), linewidth = 1.2) +
  geom_pointrange(data = dme2.2, aes(x = reactive, y = estimate, ymin = conf.low, ymax = conf.high), linewidth = 0.5) +
  ylim(-0.5, 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of Prior PI on \n Party Institutionalization") +
  xlab("Reactive ASP") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold"))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_d_4_2.jpg", dpi = 500, width = 5, height = 5) 

#reactive*pi4####
lm1 <- lm(v2xps_party ~ gwf_prior + reactive*pi4 + prior_pi + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = df)
me1 <- plot_slopes(lm1, variables = "reactive", condition = "pi4", vcov = ~ccode) 
dme1 <- me1$plot_env$dat

interplot(m = lm1, var1 = "reactive", var2 = "pi4", hist = TRUE, adjCI = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of a Reactive ASP on \n Party Institutionalization") +
  xlab("ASP PI") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold")) +  
  ylim(-0.5, 0.5) +
  geom_ribbon(data = dme1, aes(x = pi4, y = estimate, ymin = conf.low, ymax = conf.high), alpha = 0.1) 

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_d_5_1.jpg", dpi = 500, width = 5, height = 5) 

lm2 <- lm(v2xps_party ~ gwf_prior + reactive*pi4 + prior_pi + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = df)
me2.1 <- plot_slopes(lm2, variables = "pi4", condition = "reactive") 
dme2.1 <- me2.1$plot_env$dat

me2.2 <- plot_slopes(lm2, variables = "pi4", condition = "reactive", vcov = ~ccode) 
dme2.2 <- me2.2$plot_env$dat

ggplot() +
  geom_pointrange(data = dme2.1, aes(x = reactive, y = estimate, ymin = conf.low, ymax = conf.high), linewidth = 1.2) +
  geom_pointrange(data = dme2.2, aes(x = reactive, y = estimate, ymin = conf.low, ymax = conf.high), linewidth = 0.5) +
  ylim(-0.55, 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of ASP PI on \n Party Institutionalization") +
  xlab("Reactive ASP") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold"))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_d_5_2.jpg", dpi = 500, width = 5, height = 5) 
